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EXPLORING LASER-DRIVEN QUANTUM PHENOMENA FROM A 
TIME-FREQUENCY ANALYSIS PERSPECTIVE: A COMPREHENSIVE STUDY 

YAE-LIN SHEU, HAU-TIENG WU, AND LIANG-YAN HSU 


Abstract. Time-frequency (TF) analysis is a powerful tool for exploring ultrafast dynamics in atoms and 
molecules. While some TF methods have demonstrated their usefulness and potential in several of quantum 
systems, a systematic comparison among these methods is still lacking. To this end, we compare a series of 
classical and contemporary TF methods by taking hydrogen atom in a strong laser field as a benchmark. 
In addition, several TF methods such as Cohen class distribution other than the Wigner-Ville distribution, 
reassignment methods, and the empirical mode decomposition method are first introduced to exploration of 
ultrafast dynamics. Among these TF methods, the synchrosqueezing transform successfully illustrates the 
physical mechanisms in the multiphoton ionization regime and in the tunneling ionization regime. Further¬ 
more, an empirical procedure to analyze an unknown complicated quantum system is provided, indicating 
the versatility of TF analysis as a new viable venue for exploring quantum dynamics. 


1. Introduction 

Time-dependent quantum mechanics is a fundamental topic in physics, chemistry, and engineering. Tra¬ 
ditionally, dynamics of a quantum system such as lifetime and energy difference between states can be 
revealed by spectral lines in the frequency domain by using the Fourier analysis. However, the Fourier trans¬ 
form presents limited chronological information of a dynamical system. To explore chronological information 
in a quantum dynamical system, time-frequency (TF) methods are necessary. 

In recent experiment advances, TF methods such as wavelet transform combined with ultrafast spec¬ 
troscopy techniques were used to probe dynamics in molecules, solids and liquids particular 

in light harvesting complexes[8j[9, IT. In attosecond physics, the short-time Fourier transform (STFT), the 
continuous wavelet transform (CWT), Wigner-Ville distribution (WVD) (which belongs to the Cohen class 
distribution), were employed to uncover dynamical mechanisms including high-order harmonic generation 
(HHG) [HJ [T2j [13] HH H5j [16] |T7] . Lately, Hilbert transform and the synchrosqueezing transform (SST) were 
applied to investigate multiple scattering m and the dynamic origin of near- and below-threshold harmonic 
generation, respectively [19j [20] . While these studies show the usefulness and potential of each individual TF 
analysis for probing dynamics of a quantum system, several crucial issues remain when it comes to analysis 
of a quantum system with unknown complicated dynamics. First, different types of TF methods and the 
choice of window functions may provide extremely different TF representation, leading to conflicting physical 
interpretation. Second, it is difficult to select a particular TF method based on previous studies, in which 
the physical models are different, e.g ., ID [I2j [18] , 2D [12] and 3D atoms [EJ [I5j [I6j [I9j [20] , and solved by 
different numerical schemes [Hum la Eg. Clearly, TF representations derived from different TF methods 
based on different physical models cannot provide an impartial comparison. Third, in the past decades, 
several modern TF methods, e.g., the Cohen class distribution, empirical mode decomposition with Hilbert 
spectrum (EMD-HS), reassignment methods (RM), and SST have been proposed and successfully applied to 
classical macroscopic dynamical systems including molecular dynamics [21], cardiopulmonary coupling phe¬ 
nomena [22] . chronotaxic systems [23], lamb wave propagation [2T and seismic data [25] . However, Cohen 
class distribution other than the WVD, EMD-HS, RM, and different forms of SST have not been discussed 
in quantum dynamical systems. As a consequence, a comprehensive TF study in the same benchmark is 
essential. 

To this end, we take 3D hydrogen atom as a benchmark because hydrogen is one of the most representative 
systems in quantum mechanics. All simulations are based on the same numerical method (time-dependent 
generalized pseudospectral method), and different contemporary TF analysis, including the STFT, CWT, 
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Cohen class distributions, RM, SSTs and EMD-HS, are applied to study the simulated signals. The system¬ 
atic comparison of the TF methods enables to interpret physical processes in a quantum system. 

This article is organized as follows. In Section [2] we summarize the physical model of the hydrogen system 
in a strong laser field and the numerical simulation. In Section [3J we provide an overview of the TF methods, 
as well as a discussion of their pros and cons. In order to show that quantum dynamical phenomena can be 
depicted by the TF methods, we perform simulations with two different sets of laser parameters, including 
hydrogen in the regime of multiphoton ionization and tunneling ionization. The results of TF representations 
by different TF methods and a discussion of their limitations and potential applications are given in Section[4j 
In Section [5] we conclude the paper and provide an empirical procedure to analyze an unknown complicated 
quantum system. 


2. Theoretical modeling 


We simulate hydrogen dynamics in a linear polarized laser field in the framework of the electric dipole 
approximation and non-relativistic quantum mechanics. Note that for the wavelength range of approximately 
500 to 1200 nm, the electric dipole approximation is valid [2b only when the laser intensity Iq is smaller 
than 10 15 ^ 10 16 W/cm 2 . The time-dependent Schrodinger equation for atomic hydrogen interacting with a 
linearly polarized field along the z-axis in atomic units can be expressed as 


(i) 


dt 




VA ,t), 


where t/>(r, t) is the wave function at position r = (x, y , z) and at time t, z = r cos(0), and E(t) is the external 
laser field. 

To numerically solve this equation, we adopt a time-dependent generalized pseudospectral method [27j [28] 
which consists of two essential steps: (1) The spatial coordinates are optimally discretized in a nonuniform 
fashion by means of the generalized pseudospectral technique: the grid is denser near the origin and sparser 
away from the origin; (2) A second order split operator technique in the energy representation, which allows 
the explicit limitation of undesirable fast-oscillating high energy components, is used to obtain an efficient 
and accurate time propagation of the wave function. 

According to previous studies mum, dynamical phenomena such as the HHG is associated with the 
electric dipoles, i.e., the laser-driven electron oscillating around the stationary nucleus, which could be 
expressed as, respectively, the time-dependent induced dipole in the length form, denoted as ^(t), and in 
the acceleration form, denoted as d/\ (t) [28] : 

(2) d L (t) = J 'ip*(r,t)z'ip(r,t)dr. 

(3) d A (t) = J dr. 


3. Time-frequency Analysis 

Suppose the signal x(t) is composed of finite oscillatory components, that is, x(t) = J2i=i a i(t) cos(27r</>/(£)), 
and each component has a time-varying amplitude modulation (AM), ai(t) > 0, and time-varying instanta¬ 
neous frequency (IF) </>J(£), where ' stands for the first derivative, then the TF representation reflects the IF, 
which is a generalization of the notation frequency, and the localized phase information could be extracted. 
We call a,i(t) cos(27T0/(t)) an intrinsic mode type (IMT) function. More details about this kind of function 
could be found in Appendix A. 

To disclose the time-varying nature of this kind of signal, in particular the IF and AM, the TF analysis 
is a powerful approach. In the following subsections we provide an overview on classical and contemporary 
TF methods. 
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3.1. Linear TF methods. The global analysis nature of the Fourier transform is responsible for its lim¬ 
itation in extracting time-varying dynamics inside an oscillatory signal. An intuitive way to resolve this 
limitation is analyzing the signal locally; that is, we could crop a segment of finite length and apply the 
Fourier transform, and expect to observe how a signal oscillates locally. This idea leads to the STFT. We 
introduce a window function g to crop the signal at different time and perform the Fourier analysis. This 
analysis results in a time-frequency (TF) representation of the signal [51] : 

/ oo 

x(C)g*( C - 

-oo 

We call STFT 9 x (t,uj) a TF representation of the given function x and |STFT^(£, uj)\ 2 the spectrogram of x. 
Here the STFT by Eq. © is different from the conventional definition by an additional modulation factor 

e~ iujt . 

Nevertheless, temporal and frequency resolutions cannot be achieved simultaneously by the STFT, accord¬ 
ing to the Heisenberg uncertainty principal m- In other words, a wide window provides a good frequency 
estimation in the STFT at the cost of poor temporal resolution, while a narrow window has the opposite 

trade-off. In this research we follow the tradition and choose the Gaussian function with a standard deviation 

_ C 2 

of (j (the full width at half maximum of the Gaussian is 2\/21n2a) as the window function, i.e., g(Cf) = 
which leads to the Gabor transform (GT): 

_ /■£ _£\ 2 

x(Qe~~^ e -»w(c-t)d£_ 

-OO 

Here the GT we apply here differs from the usual one by an additional modulation factor e luJt . 

Based on the same local analysis idea as that of the STFT, we could analyze the momentary behavior of 
the signal x{t) by the CWT: 

(6) CWTUt,a) = J™ x(0^g* 

where g is the chosen mother wavelet and the scale parameter a controls the dilation of the window, and * 
denote the complex conjugate mm- Due to the dilation nature of the transform, the TF representation 
analyzed by the CWT has a good time resolution and poor frequency resolution at the high frequency region 
and a good frequency resolution and poor time resolution at the low frequency region. One of the commonly 
applied mother wavelet is the Morlet wavelet, and the CWT based on the Morlet wavelet is called the Morlet 
wavelet transform (MWT): 

/ oo 

x(OV^g*HC-mc, 

-oo 

-c 2 

where uj > 0, g(Q = is the Morlet wavelet and r > 0. The standard deviation for the dilated 

Morlet wavelet is a = t/uj. In Eq. (|7|), we follow the convention and use cu, which is the inverse of the 
scale parameter a in Eq. ©• Note the fundamental difference between Eq. 0 and Eq. 0 - in the GT, the 
Gaussian function is of a fixed width but in the MWT the width varies. In the MWT the window width 
a varies with uj such that the quality factor r (the inverse of the relative bandwidth) is invariant on the 
TF plane [31]. In other words, the window width a becomes smaller as uj increases, and vice versa. Since 
the frequency resolution depends on the scale, we say that the TF representation by CWT has an adaptive 
frequency resolution. 

Both the STFT and the CWT belong to the linear type TF analysis , in which the signal is characterized 
by their inner products with a preassigned family of templates with free parameters. Clearly, these TF 
representations depend on the chosen window, which might cause artificial patterns on the analysis result. 
Further, while in general the underlying structure of the signal under analysis is not known, there is no 
systematic way to design the window which faithfully reflects the structure. The above facts render the 
linear type TF analysis non-adaptive to the signal under analysis. 
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3.2. Quadratic TF methods. To resolve the non-adaptivity issue of the linear type TF analysis as well 
as finding the higher order structure, the Wigner-Ville distribution (WVD) 31, was proposed. The WVD is 
based on the concept of the autocorrelation function and is defined as 

/ oo 

x(t + C/2)x*(t-C/2)e-^dC. 

-OO 


Note that by a slight change of variable, we could view the WVD as a variant STFT which is free of the 
window choice issue. In this sense, WVD is adaptive to the signal under analysis. Also note that by a direct 
derivation, the spectrogram could be understood as follows (31]: 


(9) 


|STFT?(*, W )| 2 



WVD X (C, r?)WVD ff (C -t,n- uj)dCdrj. 


The WVD has several good mathematical properties, e.g., the signal energy is preserved in the WVD, the 
WVD is a real-valued function on the TF-plane and it provides a precise information about the chirp signal, 
like x(t) = e *M/3o+/3it+/3 2 t 2 )^ w p ere /3 0 G M, /3i > 0 and /? 2 G M. However, choosing the signal itself as the 
window function causes specific artifacts depending on the signal type. For example, when there are more 
than one oscillatory component in the signal, the interference patterns in the TF representation is inevitable, 
which might lead to mis-interpretation of the signal. 

By construction, the WVD is quadratic in the signal which could be viewed as an energy distribution. 
A direct generalization of the WVD based on imposing some constraints on the covariance structure [31] of 
the signal leads to the Cohen’s class. In other words, the Cohen’s class , which comprises all bilinear TF 
representations that are covariant to shifts in both time and frequency m and has the following form: 

(10) c x (t,u>) = JJJ x( S + C/2)x*( S -C/2)e- i ^e- i « s - f )/(^C)d^dsdC, 


where /(£, () is an arbitrary parameter function. Different parameter functions lead to different TF repre¬ 
sentations. Note that the WVD is a member of the Cohen’s class when /(£,£) = 1. Different parameter 
functions lead to different TF techniques. 

One particular technique is eliminating the interferences in the WVD when there are more than one oscil¬ 
latory component. To be more specific, we can employ a separable parameter function /(£, () = G(£)h(Q = 
Wft(£,C) in Eq. (10), where suitably chosen G(£) and h(() permit a continuous and independent control of 
the interferences in time and frequency, respectively. The corresponding representation is called the smoothed 
pseudo-Wigner-Ville distribution (SPWVD) (31] and could be rewritten as 


(11) 

( 12 ) 


SPWVD? 


/ OO POO 

g(t — s) / h(()x(s (/2)x^(s — (/2)e~ 2UJ ^d(ds 

-oo J — oo 


//: 


w\m x (s,ow h (s-t^-uj)dsd^ 


where g(Q is the inverse Fourier transform of G(£). Particularly, G(£) and h(() are even functions with 
g(0) = 1 and H( 0) = 1, where H(£) is the Fourier transform of h((). 

We mention that we could consider a variation of the Cohen’s class to get the family of the affine class 
ED, e.g., the affine WVD or the affine SPWVD. The transform considered in the affine class is the time-scale 
covariant and the TF representation also has the adaptive resolution feature. To simplify the discussion, we 
do not study the affine class transforms, but we mention that their performance is similar to those in the 
Cohen’s Class. The above-mentioned transforms are overall called the quadratic TF analysis. A higher-order 
generalization is also possible, and we refer the reader of interest to, for example, [33 . [34]. 


3.3. Reassignment Methods. As features of a TF representation computed by a linear type of TF method 
are smeared by the introduced window functions and those by the quadratic type ones are obscured by 
interferences, the RM was proposed to sharpen the resolution of the TF representation. Generally speaking, 
the coefficients of the TF representation at (t,u) is reallocated to a different point (£, u) according to a 
predefined reallocation rule [35] 136] 137] . A common choice of the reallocation rule is to assign values of a TF 
representation to the local centroids. Note that this is different from the averaging idea behind the STFT 
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(Eq. (|9|). In this study we apply the RM to the STFT and SPWVD, and called the methods RM-STFT 
and RM-SPWVD, respectively. 

The reallocation rule for the RM-STFT is derived by estimating the local centers of gravity, denoted as: 


(13) 




jR {STFT^f (t, ??)(STFT®(t, ??))*} 
|STFT®(i, 77 ) | 2 


, when STFT 9 x (t,r]) ^0 
, when STFT x (t,rj) = 0. 


(14) 


' -3{STFT2®M(STFT®M)*} 
| STFT® (*,»,) 1 2 “ 

OO 


, when STFT 9 x (t,rj) ^0 
, when STFT^(£, 77 ) = 0. 


The notations tg and dg stand for the window function tg(t ) and the first derivative of g(t) in the STFT, 
respectively. 

The RM-STFT is thus defined as 


(15) 


RM-STFT® (t, u) 



|STFT^(t, u)\ 2 5(i x (t, uj) — t, u) x (t, uj) — uj)dtduj. 


Essentially, this formula reassign the energetic contents of the spectrogram to the new location (t x ,LJ x ). As 
a consequence, the RM leads to a substantially improved resolution in the TF representation. Note that the 
reassignment rules Eq. ( [T3| ) and Eq. (14) can lead to the group delay and the IF of the bandpass filtered 
signal y(t) = STFTjjht.iod [37|. using only the unwrapped phase of the STFT^(£, co>). However these physical 
meaningful expressions are numerically inefficient. Other forms of reassignment operators can be found in 


Next we consider SPWVD. Although the SPWVD can smooth out the interferences in the WVD, the 
smoothing functions introduce artificial broadening. Applying the reassignment technique on the TF repre¬ 
sentation of the SPWVD can reduce such artifacts. Similar to Eq. (15), the reassigned representation of the 
SPWVD is 


(16) 


RM-SPWVD x (t,u) 


//: 


SPWVD^(£, oj)5(t x (t , uo) — t, cd x (t, uo) — Cj)dtduo. 


where the reassignment rule (t x ,Ld x ) is determined by the concept of expectation: 


(17) 

(18) 



sWVD x (s, £)W h(s -t,£- w)dsd£ 
?WVD x (s, - t,£ - w)dsde 


Despite RMs are intuitive techniques to sharpen the linear type and quadratic type TF representations, 
inverse routines as well as mode reconstruction are not available. 


3.4. Synchrosqueezing transform. In this section we describe the SST, which is a special RM aiming 
to address the intrinsic blurring issue in the linear TF methods. To be more precise, the SST manifests IF 
characteristic according to the reallocation rule that consists of solely the frequency information, rather than 
the centroid of the TF representation [38j 39 . In addition to sharpening the TF representation, the causality 
property of the signal is preserved in the SST, which allows the decomposition of oscillatory components 
when the signal is composed of several oscillatory components. We refer the reader who has interest in 
the theoretical analysis and other details, like reconstructing each oscillatory component and robust to 
noise, of SST to [38| 1391140] . In this section we describe the synchrosqueezed STFT (SST-STFT) and the 
synchrosqueezed CWT (SST-CWT). 
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3.4.1. SST-STFT. Consider a STFT^(£,u;) of x(t) with a window function g(t) such that supp(G(u;)) C 
[— |], where G{uj) is the Fourier transform of g and d is the smallest gap of IFs of any two consecutive 

oscillatory components. The SST-STFT is given as 


(19) 


SST-STFT 9 x (t,u> 


)-f 


STFT %{t,rj)- 


\u-* x (t, V )\ 2 


aWTT 




where & x (t,r]) is the reallocation rule and a > 0 is a controllable smoothing parameter for the resolution, 
which in practice is chosen to be small. The reallocation rule given by the following equation utilizes the 
phase information hidden inside the smeared TF representation: 


( 20 ) 


& x (t,ri) 


-idt STFTgM) 
STFT 

oo 


when STFT 9 x (t,r]) ± 0 
when STFT 9 x {t,rj) = 0. 


Note that only the frequency reassignment operator is considered in the SST-STFT, so that the causality 
property of the signal can be preserved. A slight modification of Cj{t, rf) based on the second-order information 
of the IF [41] might further improve the TF representation. Indeed, when the interested oscillatory component 
could be approximated by a chirp function mm, although the TF representation is sharpened by the SST- 
STFT, a mild spreading is inevitable. We call this mild spreading the diffusive pattern. To cope with the 
diffusive pattern, we could consider the following reassignment rule: 


( 21 ) 


f u x (t,rj) + c(t,rj)(t - t x (t,rj)) when d u i x (t,rj) ^ 0 
I (2) x (t, rf) otherwise, 


where t x (t,r ]) is defined as 

(22) i(t,r?)=t + i ^^ ( ^ ) and c(t, rj) = • 

STFT y x (t,T]) d u t x (t,T)) 

The reconstruction formulas [39] for each IMT function from SST®’ STFT (f, uj) is 

r^k C0+§ 


(23) 


Xk(t) = 3? 


Wf 

y J u k 


Ju k (t )-| 

where e<l, !R denotes taking the real part, and C g = g( 0). 


SST-STFT^ (£, uj)du 


3.4.2. SST-CWT. The definition of the SST-CWT is similar to that of the SST-STFT. For a CWT with a 
mother wavelet g such that supp(G(cj)) C [1 — A, 1 + A], with A < the SST-CWT is given by 

1 _ \u-U} x (t,r ])\ 2 

— -j=e o dr?, 

«V 7r 


(24) 


/ OO 

t ) - 3/2 CWT 9 I («,t ) ) 

-oo 


where a > 0 is a smoothing parameter and Co x (t, rf) is the reallocation rule defined by 


(25) 


Vx(t,ri) 


-idtCWT^v) 

CWT l&ri) 


when CWT^(t,? 7 ) 0 

when CWT^(£,? 7 ) = 0. 


When describing an oscillatory component with a fast varying IF, the TF representation in SST-CWT 
may show a slight spreading. While the second-order information of the IF has been discussed in the 
SST-STFT [41], it is not yet been studied. 

The reconstruction formula [39, [40] for each IMT function Xk(t) under SST-CWT is 

f pWk(t) + § 

(26) x k (t) = 5ft < Rf 1 / a -3 / 2 SST-CWT^(£, a)da 

l J"k(t)~ | 


where e«l, and R g := ^^-dg. 

The results of the SST are “adaptive” to the signal in the sense that the error in the estimation depends 
only on the first three moments of the mother wavelet instead of the profile of the mother wavelet [33 39, 
[40] . In other words, the influence of the chosen window on the associated TF representation is minimized 
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compared with the linear TF methods. In addition to the above properties, the SST is robust to several 
different kinds of noise, which might be slightly non-stationary m • An important property shared by the 
SST-STFT, SST-CWT and RM-STFT is that by taking a short window, the fast varying IF could be well 
captured. 

At the first glance, the results of the reassignment technique and SST seem to break the well-known 
Heisenberg uncertain principle, which says that the temporal and frequency resolution cannot be achieved 
simultaneously. However, in the reassignment technique and SST, we have shown that we could obtain a TF 
representation with an almost perfect time and frequency resolution. The main reason is that we focus on 
the oscillatory signals, in particular the adaptive harmonic model but not the whole L 2 space. We mention 
that the behavior of the RM and SST on the general L 2 function is an open question. 


3.5. EMD. One popular way to define the IF of a given signal x(t) is via finding the analytic representation 
of x(t) with the aid of the Hilbert transform: 

(27) x(t) = x{t) + iH(x(t)) = a(t)e l ^^\ 

where T-L(x(t)) denotes the Hilbert transform of x(t). In this equation, a{t) and <F(£) are the modulus and 
phase of x(t), respectively [42] 03]. The IF of x(t) is thus defined as the rate of the varying phase: 

d<F(t) 


(28) 


u>(t) = 


1 


27r d t 

However, the IF provided by Eq. ( [28] ) is not always meaningful. To deal with this problem, the empirical mode 
decomposition algorithm (EMD) [44] was proposed to decompose x(t) into a series of functions called the 
intrinsic mode function (IMF), x = XlfLi x k{t)> on which the Hilbert transform can be applied subsequently, 
via an algorithm called the sifting process. IMFs must satisfy the following two conditions: (i) In the whole 
data set of a signal, the number of extrema and the number of zero crossings must either be equal or differ at 
most by one; (narrow band) (ii) At any point, the mean value of the envelope defined by the local maxima 
and the envelope defined by the local minima is zero (adoption of local properties). Note that in general an 
IMF is different from an IMT function. Details of the EMD algorithm could be found in, for example, [44] . 


For each IMF Xk(t), the corresponding IF and AM amplitudes can be estimated by Eq. (27) via the 
Hilbert transform and Eq. ( |28| ). By assigning the IF and amplitudes of all IMFs on a TF plane, we obtain 
the Hilbert spectrum (HS) HS x (t,cj) (HS), which is a TF representation for x(t) determined by 


(29) 


HS x {t,u) = Ya k mu ~ Uk(t)), 


where S(cj) denotes the Dirac delta measure. Although the EMD along with HS has been applied to several 
fields, due to a number of heuristic and ad hoc elements in the EMD algorithm, it is difficult to analyze its 
accuracy and limitation. In addition, despite the solid mathematical support, the Hilbert transform might 
be limited when applied to analyze the momentary dynamics of an oscillatory signal. First, due to the slow 
decay nature of the kernel in the Hilbert transform, keeping the causality of the signal structure might be 
difficult. Second, when the signal has time-varying amplitude and frequency, in general there is no guarantee 


to get the correct analytic signal from a given real oscillatory signal by the Hilbert transform in Eq. (27). 


We mention that the part of using the Hilbert transform to obtain the AM and IF of each IMF could be 
replaced by using the SST-STFT or SST-CWT. In this case, the method is regarded as the EMD-SST-STFT 
and EMD-SST-CWT. 

A comparison of different TF methods discussed in this study is summarized in Table 1. 


4. Results and Discussions 

In this section we employ the aforementioned TF methods on the time-dependent dipole computed by the 
time-dependent generalized pseudospectral method at the ab initio level. In a strong laser field, a variety of 
dynamical processes can occur in atoms and molecules, such as ponderomotive effect, the AC Stark effect, 
HHG [46], multiphoton ionization[47] [48] 09], above threshold ionization, tunneling ionization, rescattering 
of electron wavepacket [50] EE], etc. Up to date, the STFT [H] 02] 03], CWT 04] 06] 05], WVD 07], 
and Hilbert transform 08j have been engaged in investigation of estimation of emission time and multiple 
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Table 1 Summary of TF methods in This Study 


Type 

Multiresolution 

Choice of 
Parameters 

Inverse 

Transform 

Artifacts in TF Representation 

STFT 

No 

Yes 

Yes 

broadening by the window function 

CWT 

Yes 

Yes 

Yes 

broadening by the window function 

WVD 

No 

No 

Yes 

interference patterns 

SPWVD 

No 

Yes 

Yes 

broadening by the filter function 

RM-STFT 

No 

Yes 

No 

causality is not preserved 

RM-SPWVD 

No 

Yes 

No 

causality is not preserved 

SST-STFT 

No 

Yes 

Yes 

Diffusive pattern in the fast-varying IF 

2nd order SST-STFT 

No 

Yes 

Yes 


SST-CWT 

Yes 

Yes 

Yes 

Diffusive pattern in the fast-varying IF 

EMD-HS 

Yes? f 

Yes £ 

Yes/No o 

mode mixing, etc.* 

EMD-SST-STFT 

Yes? f 

Yes J 

Yes/No o 

mode mixing 


f: The multiresolution-like behavior of the sifting process was studied in [45] . p The stopping criteria of 
the sifting process depends on tuning several parameters, o: The inversion from the HS is not possible 
for the definition (29). It is possible if we define the HS as — cjk(t)). *: It is sensitive 

to noise; a shortly existing oscillatory component destroys the whole analysis. 


scattering. However, to the best of our understanding there is no comprehensive study comparing these TF 
methods in a quantum dynamical system. 

By taking the well-studied atomic hydrogen as a benchmark, here we focus on the laser-driven hydrogen in 
the multiphoton ionization regime and in the tunneling ionization regime. When the laser photon energy Huj is 
much smaller than the ionization potential I p , the multiphoton ionization and the tunneling ionization can be 
described by a dimensionless Keldysh parameter = y/I p /2U p, where Up = Eq/ 4cJq is the ponderomotive 

potential OSH]. Similar physical dynamics can be achieved for any atom-field interaction given a fixed 
[52] . Generally speaking, multiphoton ionization of atoms become dominant when y^ 1, while tunneling 
ionization is predominant when y k 1 [29]. In the following we perform simulations with laser parameters 
in these two regimes. 

4.1. Multiphoton Ionization Regime. In the first simulation, the laser field parameters are arranged such 
that the Keldysh parameter is y k = 3.07, indicating that multiphoton ionization is the major mechanism. 
The laser wavelength is 880 nm, which corresponds to ujq ss 0.05178273 in the atomic unit (a.u.), and the 
laser intensity of Jo = 10 13 W/cm 2 , which corresponds to a field amplitude Eq ss 0.0169 (a.u.). Note that 
for this y k value, both tunneling and multiphoton ionization can occur, but the latter is dominant. The 
laser field (Fig. 0 a)) E(t) = EoF(t) sin(co’o^) has a profile of F(t) = sin 2 (7r£/(nT)), where n = 60 is the pulse 
length measured in the optical cycle (T = 2tt/ujo). 

The computed induced dipole in the length form dp{t) is presented in Fig. [ljb). The power spectrum 
[28] . computed by the squared Fourier spectrum, of the laser field and dp{t ) are shown in Fig. [lj c )and[ljd), 
respectively. Both the length and the acceleration forms of dipole moment present the same detail structures 
in their power spectra [28] . Here we present the analysis using dz,(t) for the multiphoton ionization case. 

Although the profiles of dp(t) and the applied laser field look similar in the time domain, they are different 
in the frequency domain. For example, while the power spectrum of the laser field has only one peak located 
at cjo, that of dp(t) reveals odd harmonics due to the parity symmetry [29]. However, the meaning of the 
substructures within the odd harmonics and their corresponding dynamics are unclear. 

To unveil the dynamics of dx,(t), first we apply the conventional linear and quadratic TF methods, as 
shown in Fig. [2] The TF representation of the GT, computed by Eq. (|5|) with a = 57.94 a.u., and the MWT, 
computed by Eq. <§ with r = 30, are displayed in Fig. [2|a) and Fig7|2|b), respectively. Both the GT and 
MWT depict separate broad lines regarding the odd harmonics in the HHG process. While frequencies of 
the extracted components are consistent with the information shown in the Fourier spectrum in Fig. [ljd), 
the GT and MWT further capture the momentary behavior of each component. For a fixed resolution in the 
GT, low frequency components could be lost given an insufficient window width. In the TF representation 
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(C) (d) 

Figure 1. The simulation of laser-driven hydrogen in the multiphoton ionization regime. 
The laser wavelength is 880 nm, and the laser intensity of Jo = 10 13 W/cm 2 , corresponding 
to the Keldysh parameter of 7 k = 3.07. (a) The laser profile, (b) The induced dipole 
moment dz,(£). (c) The power spectrum of the laser profile, (d) The power spectrum of 
d,L(t), respectively. Note that the laser profile and di,(t) are very similar, yet very different 
in their spectral components. 


of MWT, the same components are expected to be captured but the representation should be different 
from the GT due to the dilation nature of the transform; that is, the adaptive frequency resolution of the 
transform. Indeed, the frequency resolution below the 7th harmonic is improved compared with Fig. |2|b), 
while harmonics on the upper TF plane remain broaden. In addition, the frequency-dependent weighting 
in y/uj enhances the representation for high order harmonics. Note the width of the window function g(t) 
is uniform in the GT, and increases as uo increases in the MWT, as shown by the vertical lines between 
neighboring harmonics in the upper TF plane. The crescent shape distribution near the boundaries of the 
TF plane is caused by the boundary effect in CWT - for a low frequency, the corresponding scale is large 
and the influence of the window cut-off near the boundary becomes apparent. The boundary effect is also 
inevitable in the GT, like the vertical artificial pattern in the very beginning and end. However it is less 
dominant. 

As discussed above, while the window is avoided in the WVD, a strong interference pattern is inevitable 
[51] . The TF representation of the WVD, shown in Fig.[2jc), reveals components between the odd harmonics, 
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which is inconsistent with the features in the Fourier spectrum and violates the parity symmetry in physics. 
Moreover, harmonics above 15, of which intensities are weak, are not observed within the range of colorbar. 
Interference along the time-direction is also observed. As there are more oscillatory components and more 
overlaps caused by x(s + £/2)x*(s — C/2) in Eq. (J8| as time increased, interference becomes stronger, as 
shown on the right-half TF plane in Fig. [2jc). Despite generating inaccurate information, note that in this 
simulation the WVD suggest the features of the TF representation, i.e., slow-varying components, regardless 
of a window parameter. This suggest that the window width in the linear type TF methods should be large 
enough such that patterns depict the harmonics as that performed by the WVD. 


In order to remove the interference, we employ the SPWVD with the filtering functions g(s) = e 2cj a 2 with 

-g 2 

G g = 100 and H(£) = e 2cr H 2 with gh = 0.025. Note that g( 0) = 1 and H( 0) = 1. The filtered result is given 
in Fig. [ljd). Although the even harmonics are eliminated, temporal interference cannot be fully removed. 
The introduced filters also worsen both temporal and frequency resolutions, as well as generating unexpected 
features. Note that the dynamic range (colorbar) is increased to reveal weak details. 

By reassigning the energy distribution of the linear type TF methods and the SPWVD to their local 
centroids, the TF representations are sharpened. The results of reassigned GT (RM-GT) and reassigned 
SPWVD (RM-SPWVD) are shown in Fig. |3ja) and (b). In Fig. [3ja), the TF representation of the RM-GT 
significantly improved the resolution, revealing clear and distinct odd harmonics, as well as their subtle 
variation. Note that it is the square of the GT representation that is reassigned, therefore the range of the 
color bar differs from Fig. la). 

In Fig. [3jb), although the broadening that comes with filtering functions are removed and the TF repre¬ 
sentation depicts odd harmonics similar to those in Fig. |3ja), the interference artifact cannot be eradicated. 
We further employ the SST on the GT (SST-GT) and MWT (SST-MWT), in which the allocation rule is a 
good approximation of the IF of each IMT. The results shown in Fig.[3^c) and (d) also depict similar features 
as that in the RM-GT. Note that the SST-MWT preserves the adaptive resolution, weighted frequencies, 
and the boundary effect of the CWT. 

In Fig. § a), (c) and (d), we show that the IMT functions can be decomposed without ambiguity by 
the RM-GT, SST-GT, and SST-MWT. Although the RM can illustrate the distinct odd harmonics and 
frequency shift of the AC Stark effect by finding the local centroids, there is no routine for reconstructing 
the harmonics. SST, on the other hand, comes equipped with an inverse method (See Eq. (23) and Eq. (26)). 

To explore the physical meaning of the shifting in Fig. |4'a)-(c), we analyze it using the Floquet method 
[53] , which has been extensively used in chemical physics [54] [55]. Details of the Floquet method can be 
found in Appendix B. In Fig. [4jd), the blue, green, red, and cyan lines denote the energy difference of ls-2s, 
ls-2p x , ls-2p z , and ls-2p y , respectively, computed by the Floquet method, in the unit of ujq. Due to the 
breaking of the spherical symmetry by the laser field, the energy levels of 2s, 2p x , 2p y , and 2p z shift and split. 
The quasi-energies derived by the Floquet computation and the curve by the SST-GT are superimposed for 
comparison. We found that the spectral line on the TF representations obtained by the SSTs not only 
demonstrate the AC Stark effect but also depict the selection rule, i.e., only the ls-2p z energy difference is 
presence. 

It is worth mentioning that figures in Fig. [3] illustrate a frequency shift at the beginning cycles of the 
7th harmonic, which could be the AC Stark effect. In the following context, we discuss about whether the 
IF components obtained by the reassigned values have a physical meaning, or they are simply clamped by 
artificial processing. 

Details around the 7th harmonic are enlarged in Fig. [4] for the (a) RM-GT, (b) SST-GT and (c) SST- 
MWT. Note that the colorbar in Fig.[4^a) is different from from that in Fig. B b ) and (c), since in the RM-GT 
the squared TF representation is reassigned. Despite different intensities, the RM-GT, SST-GT and SST- 
MWT depict a similar shifting trend descending from the 7.241-th harmonic to the 7.200-th harmonic, 
which corresponds to the AC Stark effect. The 7.241-th harmonic at the beginning cycles corresponds to 
the energies for ls-2p transition (jAl — = 7.241 in a.u.), manifested as a small peak overlapping the 

H7 of the power spectrum in Fig.JlJd). The intensity of such spectral line is small because it arises from the 
near resonance absorption (not resonance absorption). 
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Figure 2. TF representations in the multiphoton ionization regime by the (a) GT, (b) 
MWT, (c) WVD, and (d) SPWVD. The lines indicating odd harmonics in the GT and MWT 
are subject to broadening issues arising from a window. Although the broadening artifact can 
be alleviated by the WVD, the evoked interference artifacts in the TF representation result 
in incomprehensible analysis. By applying additional filters, the SPWVD can moderate the 
interference patterns, at the price of broadening of the features in the TF representation. 


Fig. [4jb) and (c) suggest that the IFs shown by different SST methods are independent of the chosen 
linear TF method. Based on the independence of the chosen SST methods and the intimate matching of the 
analysis result and the theoretical prediction shown in Fig. Qd), we conclude that the decomposed IMTs 
are physically meaningful. 

The ionization process for the laser field in Fig. [T|a) is summarized in Fig. [ 5 ] In Fig.[5|a), when the laser 
field is small, the spectral line 7.241 aroused from the ls-2p near resonance absorption caused by the atomic 
structure. As the laser intensity gradually increases, the energy levels of 2s, 2p x , 2p y , and 2p z shift and split 
due to the breaking of the symmetry by the electric field, which is regarded as the AC Stark effect. When 
the laser intensity increases furthermore, the high-order harmonic process is induced and the dynamics is 
dominated by the transition between dressed states |M, N) formed by the electron state M and the photon 
state V, as illustrated in Fig. [5|b). 
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Figure 3. TF representations in the multiphoton ionization regime by the (a) RM-GT 
and (b) RM-SPWVD. The broadening caused by the window in the GT and filter functions 
in SPWVD is removed. Note that the frequency shift in the beginning 10 cycles corresponds 
to the AC Stark effect. TF representations in the multiphoton ionization regime by the (c) 
SST-GT and (d) SST-MWT. The broadening issue caused by the window in GT and MWT 
are removed. Note that the frequency shift in the beginning 10 cycles corresponds the AC 
Stark effect. 


A summary for comparison of TF methods in the multiphoton ionization regime is provided in Table 2. 
Note that in this table we mention that the odd harmonics and the AC Stark effect cannot be obtained by 
EMD with the HS. In fact, we could not yet find physical meaning for the series of IMFs decomposed by the 
EMD in the multiphoton ionization regime. 


4.2. Tunneling Ionization Regime. In the second simulation, we discuss the TF representations for 
signals in the tunneling ionization regime. The laser wavelength in this subsection is ujq ss 0.05696100 in 
atomic units (a.u.), corresponding to 800 nm, and the laser intensity is /o = 3.5 x 10 14 W/cm 2 . The Keldysh 
parameter is 7 k = 0.57, suggesting that the tunneling ionization is dominant. The laser field has 10 cycles 
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Figure 4. The AC Stark effect revealed by the (a) RM-GT, (b) SST-GT and (c) SST- 
MWT. (d) Comparison of the frequency shift caused by the AC Stark effect computed by 
the SST-GT (gray scale) and the Floquet method. 

Table 2 Comparison of TF methods in the multiphoton ionization regime 
Phenomenon Discrete Odd Harmonics The AC Stark Effect 


GT 

Yes (blurred) 

Yes (blurred) 

MWT 

Yes (blurred) 

Yes (blurred) 

WVD 

No 

No 

SPWVD 

Yes, for 

harmonics > 5 (blurred) 

Yes (blurred) 

RM-GT 

Yes 

Yes 

RM-SPWVD 

Yes 

Yes (blurred) 

SST-GT 

Yes 

Yes 

SST-MWT 

Yes 

Yes 

EMD-HS 

No 

No 


In this study, the GT and MWT are examples for the STFT and CWT, respectively. 
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Figure 5. Illustration for the (a) AC-Stark effect and (b) the high-order harmonic process. 


and its profile is ramped on according to 


(30) 


Eo(t) = 


sin2 (w) 

when 

0 < t < T 

1 

when 

T <t <9T 

sin2 (w) 

when 

9T < t < 10T, 


which leads to the laser field EoF(t)sm(uot), as shown in Fig. (6^a). We adopt the dipole moment in 
acceleration form in the case of tunneling ionization, as displayed in Fig. [6^b). The acceleration form can 
provide clearer information for the high order harmonics than the length form, because the differential 
operator in Eq. ©> acts as a high-pass filter. 

TF representations for the linear type and quadratic type of transforms are presented in Fig. [7| including 
the GT (Fig. |7|;a)), MWT (Fig. 0b)), WVD (Fig. 0c)), and SPWVD (Fig.0d)). In Fig. 0c), we observe a 
periodic repetition of arches, suggesting the chirp-like dynamics of the attosecond radiation, as predicted by 
the standard semiclassical model [29 . To present such feature in the linear type TF representations, a short 
window width and a small quality factor r are needed for the GT and the MWT, respectively. (Applicability 
of a short window width in the IMT models deserves further investigation.) Here we select a — 1 (a.u.) 
for the GT and r = 6 for the MWT. TF representations in the GT and the MWT illustrate that the arch 
structure repeats every 0.5T. The inner structure after the 2nd cycle is indistinct because of the broadening 
caused by the window. While the GT cannot resolve information below the 15th harmonic due to the small 
window width, the MWT can describe the near and the below threshold harmonics (low order harmonics) 
by its merit of adaptive resolution. WVD, on the other hand, can characterize the chirp very well without 
the need of window parameter and the broadening that comes along. However, as described in previous 
subsection, WVD introduces strong interferences in both temporal and frequency directions. Despite in 
SPWVD, where we use the parameters a g = 1 and an = 0.13, it does not seem possible to eliminate the 
interference that mingled with intrinsic structure inside the main arches. 

We further apply the RM techniques by reassigning the TF representations by the rule of centroids. In 
Fig. [IJa) and Fig. [8](b) , we show the results of the RM-GT and the RM-SPWVD, representatives for the 
linear and quadratic type methods, respectively, can address the broadening issue in the original methods 
and provide distinct inner arches as well as other substructures. Note that in Fig. |8jb), the arch in the 
duration of 2 — 2.5T and around the 45th harmonic is an example of interference that is difficult to remove. 
It is because that the RM techniques re-shuffle the original TF representation based on its local center of 
mass, and hence preserves the intrinsic artifact in the adopted transform. 

By a reallocation rule that contains phase information of a linear type TF transform, the SST can enhance 
the TF representation by linear type TF transforms in a manner similar to that by the RM technique, 
as shown in Fig. |8|c) and (d). Both the SST-GT and the SST-MWT can depict similar structure with 
appropriate the STFT and CWT parameters. Note that as mentioned previously, the SST-GT inherited the 
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(a) (b) 

Figure 6 . The simulation of laser-driven hydrogen in the tunneling ionization regime. 
The laser wavelength is 800 nm, and the laser intensity of 3.5 x 10 14 W/cm 2 , corresponding 
to the Keldysh parameter of 'Jk = 0.57. (a) The laser profile, (b) The induced dipole 
moment 


fixed resolution feature in the STFT, and the SST-MWT maintains the adaptive resolution in the CWT. By 
comparing Fig. J8^c) with Fig. [8jc), we see that the resolution in the RM-GT is more sharpen than that of 
the SST-GT. One of the reasons is that while both temporal and frequency reassignments are considered in 
the conventional reassignment method, temporal reassignment is not taken into account in the SST (both 
SST-GT and SST-MWT). In addition, the reallocation rule Eq. (20) is only the first order approximation 
for the IF, resulting in diffusive pattern particularly for fast-varying chirp signal. For the simulation in 
the multiphoton ionization regime, such diffusive pattern is negligible as the IF changes slowly. By using 
Eq. (21), a second-order reallocation rule for the SST-GT, the diffusive pattern can be further concentrated, 
as is shown in Fig. [8je). 

To shed some light on this TF representation for the tunneling ionization mechanism, we follow the 
standard semiclassical approach suggested independently by Corkum [56] and Kulander et al m- Here the 
electric-field force corresponding to the applied laser field in a.u. is F z = — E^E{t)z , where z is the unit 
vector in the z—direction. We assume the initial position and initial velocity is both zero [29] . 

The trajectory of electrons released between IT and 1.5T by using the extended semiclassical approach 
is denoted by red circles in Fig. [8^f). This result is superimposed for the sake of comparison with the TF 
representation of the SST-GT. The trajectory suggests that the largest arch in Fig. [8] is associated with the 
first return time of the electron wave packet released in between IT and 1.5T of the laser field, and the 
smaller arch between 2T and 2.5T, is the second return. 

Although the SST depicts the first and multiple returns clearly, the sequence of these returns cannot be 
known. To retrieve the first and second return quantitatively, Risoud et al propose the application of the 
Hilbert transform on the dipole moment with the harmonics lower than the ionization potential filtered PI- 
While there are more than one component simultaneously exist at each time, we consider the EMD scheme 
to decompose the filtered dipole moment. (In this case, harmonics below the 10th harmonic are filtered as 
they correspond to the bound part of the atomic spectrum.) The stopping criterion in the EMD scheme is 
that the residue becomes monotonic. The HS for all modes of IMF are displayed in Fig. [9^a). Note that a 
median filter with a window length of 51 points is applied on all modes as the implementation of Eq. (28) 
can introduce numerical errors. Despite there is no theoretical base for the IMFs by the EMD and the result 
is not consistent in the multiphoton ionization case, in this case we see that the first two modes can be 
associated with the first and second returns. Fig. j9^b) presents the HS for the first two IMFs and Fig. j9^c) 
compares the IMFs with the representation of the SST-GT. In Fig. [9](b) , the IMFs are visually emphasized 
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Figure 7. TF representations in the tunneling ionization regime by the (a) GT, (b) MWT, 
(c) WVD, and (d) SPWVD. The repeating arches suggest excitation of the harmonics is 
strongly correlated with the laser field. In the representations of the GT and MWT, the 
arches are broaden by the window, leading to the difficulty to differentiate the intricate 
structures inside the arches. The representation of the WVD is obscured by the interference 
pattern, and it is difficult to remove the interference pattern with filters in the SPWVD. 


by additional Gaussian functions in both temporal and frequency domain. Note that the IF revealed by the 
Hilbert transform is meaningful only if the IMF is consist of a single component, which is not guaranteed 
in the EMD. For example, for time greater than 2.5 cycles, the HS of the first IMF is distorted because of 
more wave packet returning from previous emissions. 

While the Hilbert transform is not theoretically suitable for this kind of oscillatory signal with time-varying 
AM and IF, we illustrate the IF of these decomposed modes by the SST-GT. We employ the SST-GT on 
the first two IMFs, as presented in Fig. |9jd) and (e). Results in Fig. [9jd) and (e) suggest that expressing 
the IF of the IMF by the SST is more stable than using the Hilbert transform. Note that no median filter 
is necessary for each IMFs in the SST-GT analysis, which reduces the possible artifacts. 

A summary for comparison of TF methods in the tunneling ionization regime is provided in Table 3. 
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Figure 8 . TF representations in the tunneling ionization regime by the (a) RM-GT 
and (b) RM-SPWVD. The broadening artifacts caused by the window in both methods 
are eliminated, but it is difficult to remove the interference pattern in the SPWVD. TF 
representations in the tunneling ionization regime by the (c) SST-GT and (d) SST-MWT. 
The broadening artifacts caused by the window in both methods are alleviated, (e) TF 
representation in the tunneling ionization regime by the second order SST-GT. Note that 
the diffusive pattern in the SST-GT is modified for the arches, (f) Comparison of the TF 
representation of the SST-GT and a trajectory of the electron released between the IT and 
1.5T (denoted by red circles) computed by the standard semiclassical approach suggested 
independently by Corkum [56] and Kulander et al. m- 





















Harmonic Order Harmonic Order 


18 


YAE-LIN SHEU, HAU-TIENG WU, AND LIANG-YAN HSU 




Figure 9. (a) TF representation in the tunneling ionization regime of the EMD-HS 

algorithm, which contains all IMF modes extracted by EMD. (b) TF representation in 
the tunneling ionization regime of the first two IMFs extracted by EMD. To enhance the 
visualization, the curves are thickened, (c) Comparison of the first two IMFs with the TF 
representation by the SST-GT. According to the semiclassical trajectory, the first IMF is 
the first return and the second IMF is the second return, (d) The first IMF and (e) the 
second IMF analyzed by the SST-GT. 
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Table 3 Comparison of TF methods in the tunneling ionization regime 


Phenomenon 

First Return 

Second Return 

Sequencing of 
Multiple Returns 

Near and Below 
Threshold Harmonics 

GT 

Yes (blurred) 

Yes 

(blurred) 

No 

No 

MWT 

Yes (blurred) 

Yes 

(blurred) 

No 

Yes (blurred) 

WVD 

No 

No 


No 

No 

SPWVD 

Yes (blurred) 

No 


No 

No 

RM-GT 

Yes 

Yes 


No 

Yes 

RM-SPWVD 

Yes 

No 


No 

No 

SST-GT 

Yes (diffused) 

Yes 


No 

No 

2nd order SST-STFT 

Yes 

Yes 


No 

No 

SST-MWT 

Yes (diffused) 

Yes 


No 

Yes 

EMD-HS 

Yes (mode mixing) 

Yes 

(mode mixing) 

Yes 

No 

EMD-SST 

Yes (mode mixing) 

Yes 

(mode mixing) 

Yes 

No 


5. Conclusions 

This paper provides a benchmark of analyzing time-dependent quantum systems by applying several con¬ 
temporary TF methods. Within the same physical model and computational scheme, different features of TF 
representations provided by different TF methods may lead to conflicting interpretation in physics. In the 
multiphoton ionization regime, linear TF methods with an appropriate window function, can vaguely depict 
discrete odd harmonics in the HHG process, yet the detail features such as the AC Stark effect cannot be 
revealed. In addition, in the tunneling ionization regime, the first-return trajectories are roughly illustrated, 
but the multiple returns are indistinct because of the broadening induced by the window. While present little 
broadening, TF representations by the WVD suffer from artificial interference patterns and lead to incorrect 
interpretation of physical mechanisms, such as even harmonics in the multiphoton ionization regime and false 
trajectories in the tunneling ionization regime. To eliminate the artifacts mentioned above, several modern 
TF analysis methods are introduced in quantum dynamics for the first time. However, Cohen class distribu¬ 
tions such as the SPWVD can only remove portion of interference patterns. Reassignment techniques and 
the SST can depict accurate dynamic process in the two regimes based on the principle of local centroid and 
instantaneous frequency, respectively. However, there is no inverse transform in the reassignment techniques. 
On the other hand, the SST preserves signal causality and allows mode reconstruction. We demonstrate that 
after removing the broadening, the AC Stark effect in the multiphoton ionization regime and the multiple 
returns in the tunneling ionization regime are revealed. Finally, we relate the IMFs from the EMD to the 
multiscattering of the electron energy distribution in the tunneling ionization regime. 

We believe that in addition to the atomic hydrogen system, the contemporary TF methods have potential 
for exploring other complicated quantum systems. In case of analyzing dynamics in an unknown quantum 
system, this paper suggests the following procedure. First, features obtained by the WVD can be used to 
determine the parameter for the linear type methods. Then the SST can be applied to obtain TF represen¬ 
tations for further interpretation. We hope that this research can serve as a cornerstone in applications of 
TF analysis for fields including but not limited to, attosecond physics, nuclear magnetic resonance, ultrafast 
dynamics in atoms and molecules. 
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Appendix A. Adaptive harmonic model 

Fix constants 0 < e « 1, 0 < d < 1, c 2 > ci > e and c 2 > cs > e. Consider the functional set Qci,c 2 ,c 3 ^ 
which consists of functions x(t) = Aft ) cos(27r 4>(t)) G C 1 (M) D L°°(M) where the following conditions hold: 

' A G C 1 (M) D L°°(R), (f>eC 3 (R), 

inf temA(t) > c lt inf te i </>'(£) > c u 

< sup teR A(t) < c 2 , sup teR </>'(£) < c 2 , sup teR \<t>"(t)\ < c 3 , 

|^4' (t) | < ecf) f (t), |</> /// (£)| < (£) for all £ G M, 

We call x = A(£) cos(27T0(£)) G Qg 1,C2,C3 an intrinsic mode type (IMT) function, </>(£) the phase function, 
the derivative of the phase function (f k {t) is called the instantaneous frequency (IF) and Aft) is called the 
amplitude modulation. When \(j)"ft)\ ecj)'(t) for all t G M, we call the signal x oscillatory with slowly 
varying IF ; otherwise we call it oscillatory with fast varying IF. Locally an IMT function with slowly varying 
IF behaves like a harmonic function and an IMT function with fast varying IF behaves like a linear chip 
function. IMT function serves as a mathematical formula for the IMF considered in the EMD algorithm, 
i.e., the IMT functions contain the properties defined in the IMF, but the reverse inclusion does not hold. 

To describe an oscillatory signal, the following adaptive harmonic model is commonly considered in 
the time-frequency analysis literature [38] (39J 00]. Fix constants 0 < e < 1, d > 0 and c 2 > c\ > 0. 
Consider the functional set Q^ 02 ’ 03 , which consists of functions in C 1 (M) nL°°(M) with the following format: 

xft) = Yl?i=i x k(t), where K is finite and x k (t) = A^(t) cos(27r<^(£)) G Q^ 1 ’ 02 ’ 03 . Let 4> r k (t) > 0 / fc _ 1 (£). In the 
STFT-SST, successive IMTs for all time t have to be separated by at least d, i.e., (f k (t) — > d. 


Appendix B. The Floquet Method 

The total Hamiltonian of atomic hydrogen can be expressed as H(t) = Hq + idi(t), where Hq is the 
unperturbed hydrogen Hamiltonian and H\(t) = — zE(t ) describes the laser-atom interaction. In the time- 
dependent generalized pseudospectral simulation, the time period T is fixed but the field amplitude varies. 
To calculate the shifting of orbital energies, we separate each cycle and take its field amplitude, and then 
assume the hydrogen in a laser field with the particular constant amplitude. For example, for the third 
cycle in Fig.[lja), the maximal amplitude is 4.86 x 10 -4 , and we assume hydrogen in the laser field with the 
constant amplitude. 

When the Hamiltonian is periodic with a time period T, i.e., H(t) — H(t + T), we can apply the Floquet 
theory [53] and derive a time-independent infinite-dimensional eigenvalue matrix equation |2I] as follows: 

(31) X! ((nlmk\H F \n , l , m , k , ))<p n J m ' k ' = e x ^ k , 

nlmk 


Here ((nlmk\HF\n'l'm'k f )) denotes the matrix elements of the time-averaged Floquet Hamiltonian over a 
period T, where the outer bra-ket notation refers to the inner product over t. We use the orbitals of an 
unperturbed hydrogen as a basis set since we only consider a weak field regime. As a result, n, /, and 
m represent the principal quantum number, angular momentum quantum number, and magnetic quantum 
number of hydrogen, respectively, k is an integer. 

The matrix element in Eq. (31) can be evaluated as 


(32) 


((nlmk | Hp \ ml I'm' k')) 

— -k^nlm^nn' &ll' &mm' &kk' 

- (nlm\z\n'l'm') + 4,/e'-1) 

T khiu5 nn ' 5n' 5 rnrn ' 5 kk ' 


The values of (nlm\z\n'l'm') can be computed analytically. For examples, in atomic units, (210|2(|l / 0 / 0 / ) = 
0.7449, (310|z|lW> = 0.2983, (210|^|3W) = 0.5418, and (310|z|2W) = 1.7695. 
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